why the additional comma if its the last element of the list?

oops :)

Not sure why one would need to init with a diagonal or spherical PLN covariance model, when the goal is to fit a PLNnetwork model, that is, one with sparse precision/covariance matrix... But if used for completeness/compatibility, that is fine.

I have exactly the same question. Diagonal covariance matrix would result in an empty graph (with only isolated nodes).

@maddyduran and I seem to remember issues when initializing PLNnetwork when you have fewer samples than species - if you use "fixed" there's a call to solve() that can throw an exception. We figured that "spherical" or diagonal would avoid this (skipping that call) and move on to estimating a penalized inverse covariance matrix even with a spherical or diagonal init.

For .logfactorial_torch(), shouldn't the final term log(pi)/2 be torch_log(pi)/2 ?

yes, though not actually sure what would happen under the hood here in terms of evaluation?

Why not use .logfactorial_torch() instead of .logfactorial() ? Because the computation has already moved back to the CPU at this point ?

Good question. Perhaps it would be better to defer that so we can use .logfactorial_torch() instead

Is there a reason to remove the convergence check on the parameter values and to keep only the one on the ELBO value ? This will speed up the algorithm (less conditions to satisfy) but may cause the nlopt and torch implementations diverge in the result they produce (not a bad thing per se, but something we need be aware of).

I defer to you on that - we have been generally happy with larger default values of ftol_rel, triggering convergence on status 3 most of the time.

I may be mistaken, but get_cov_mat() appears to be defined but never used anywhere in compute_vcov_from_resamples(). Is it necessary ?

Vestigial code from debugging - can be removed

I have exactly the same question. Diagonal covariance matrix would result in an empty graph (with only isolated nodes).

No comparison to a null Poisson model in the general post-treatment of PLN families. Is it to improve speed / because it's not required for the post-treatments ?

In our applications, the call to glm.fit inside of nullModelPoisson() can sometimes throw an exception, which ruins a perfectly good PLN fit! Since it seemed to us that nullModelPoisson() was something one only needed to do when (optionally) computing the approximate R2, we thought this call was superfluous.

I have no strong opinion on this one as @jchiquet was the one who wrote this part but is there a reason prefer lapply to future_lapply (one less dependency ? simpler to use ?). A nice thing about future_lapply is that it is backend-agnostic and can be used for several parallalelization paradigms.

Ah, yes, agree it would be wonderful to keep this. The issue is that on machines that use OpenBLAS with a multithreaded backend, using future can deadlock the session. A workaround is to wrap calls to future with something like this:

old_omp_num_threads = as.numeric(Sys.getenv("OMP_NUM_THREADS"))
if (is.na(old_omp_num_threads)){
   old_omp_num_threads = 1
}
RhpcBLASctl::omp_set_num_threads(1)

old_blas_num_threads = as.numeric(Sys.getenv("OPENBLAS_NUM_THREADS"))
if (is.na(old_omp_num_threads)){
  old_blas_num_threads = 1
}
RhpcBLASctl::blas_set_num_threads(1)

Then you do work with future and then:

RhpcBLASctl::omp_set_num_threads(old_omp_num_threads)
RhpcBLASctl::blas_set_num_threads(old_blas_num_threads)

We didn't add this because we didn't want to add a new dependency on RhpcBLASctl to the package, but you could do if you want to be able to do linear algebra inside of functions called by future

Should we add a backend element (set to the appropriate value) to config_default_nlopt and config_default_torch so that they are self-aware and that we can use
if (config$backend == "torch") {...}

That would be better. Would also be good to be able to specify the torch device (e.g. "mps", "cuda", etc)

Same as previous comment

I'm mixed about this: to avoid potential problems, it would be better to allow only covariance = "fixed" (especially since f801ca6) reverts back to PLNnetworkfit <-R6Class(inherits = PLNfit_fixedcov)

@maddyduran and I discussed this and we think the use case for covariance = "spherical" or "diagonal" in PLNnetwork is when you have fewer samples than you have species. IIRC, the issue was the call to solve() in "fixed" optimize() function.